Spin waves and local magnetizations on the Penrose tiling 
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We consider a Heisenberg antiferromagnet on the Penrose tiling, a quasiperiodic system having an 
inhomogeneous Neel-ordered ground state. Spin wave energies and wavefunctions are studied in the 
linear spin wave approximation. A linear dispersion law is found at low energies, as in other bipartite 
antiferromagnets, with an effective spin wave velocity lower than in the square lattice. Spatial 
properties of eigenmodes are characterized in several different ways. At low energies, eigenstates 
are relatively extended, and show multifractal scaling. At higher energies, states are more localized, 
and, depending on the energy, confined to sites of a specified coordination number. The ground state 
energy of this antiferromagnet, and local staggered magnetizations are calculated. Perpendicular 
space projections are presented in order to show the underlying simplicity of this "complex" ground 
state. A simple analytical model, the two-tier Heisenberg star, is presented to explain the staggered 
magnetization distribution in this antiferromagnetic system. 

PACS numbers: 71.23.Ft, 75.10.Jm, 75.10.-b 
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I. INTRODUCTION 

Theoretical investigations of the Heisenberg model in- 
volving localized spins on a quasiperiodic structure are 
motivated by experiments performed on a class of mag- 
netic quasicrystals, the icosahedral ZnMgR (R: rare- 
earth) alloys^. Short range antiferromagnetic correla- 
tions and anomalous slow dynamics have been observed 
at temperatures below 4K in ZnMgHo, and below 5.8 
K in ZnMgTb 2 . These experiments raise the question 
of what types of magnetic ordering can exist in systems 
with quasiperiodic structural order. Several theoretical 
studies of both classical^ and quantum££i£i£ models of 
interacting local moments in two dimensional quasiperi- 
odic geometries have been carried out. These calcula- 
tions show ground states characterized by a complex in- 
homogeneous magnetic order. This quasiperiodic ground 
state, and magnon modes are investigated in detail in the 
present paper. 

We consider spins located on vertices of the Penrose 
tiling, which was introduced^ as an example of a deter- 
mistic nonperiodic structure and later found to be rele- 
vant for quasicrystals^. The Penrose tiling is the planar 
version of the three dimensional icosahedral tiling used in 
the description of a large number of quasicrystalline al- 
loys. At T=0, since there is no frustration in the model, 
one expects a Neel ordered ground state, with equal and 
oppositely directed sublattice magnetizations. We will 
consider the case of S = = , the case of greatest theo- 
retical interest, as quantum effects are expected to be 
strongest. An interesting related problem concerns the 
entanglement properties of spins in such hierarchical en- 
vironments, a question that is currently receiving atten- 
tion in view of the applications to quantum computing. 

This paper follows an earlier Brief Report, where 
some early results were presented. In this paper we de- 
scribe the method used in some detail, and we present a 
more complete description of the ground state properties, 
spin wave spectrum, magnon velocity, and eigenmodes. 



The numerical calculations are done within the linearized 
spin wave (LSW) approximation, using periodic approx- 
imants of the Penrose tiling. We show as well that the 
magnetization results can be qualitatively explained in 
terms of a simple analytic model of tree clusters. 

In Sec. II we provide some background on the Heisen- 
berg antiferromagnet as well as the Penrose tiling for 
those not familiar with this quasiperiodic structure. 
Sec. Ill presents the method used, and describes the finite 
samples that we studied numerically. Sec. IV presents re- 
sults for the magnon spectrum and density of states. In 
Sec.V we give the results for the ground state energy, 
and discuss the inhomogeneous local staggered magne- 
tizations. Sec. VI presents a simple analytical model to 
explain the numerical results. Sec. VII resumes the main 
conclusions. 



II. HEISENBERG HAMILTONIAN AND THE 
PENROSE TILING 




FIG. 1: A finite patch of the perfect Penrose tiling. 
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A. The Heisenberg Hamiltonian. 

We consider the nearest-neighbor antiferromagnetic 
spin-i Heisenberg model, 

H = jJ2s i -S j , J>0, (1) 

(hi) 

where are pairs of linked vertices of the Penrose 

tiling. The antiferromagnetic coupling J is taken to be 
the same for all bonds. The tiling is bipartite, in that 
sites belong to one of two sublattices A and B, and the 
J-terms couple a spin on sublattice A to spins lying on 
sublattice B and vice-versa. This property ensures that 
there is no frustration, i.e., if one considers classical spin 
variables, the ground state is one for which all bonds are 
"satisfied" - with all the A-sublattice spins pointing in 
one direction and B-sublattice spins pointing in the oppo- 
site direction. This is our reference state with long range 
staggered magnetic order for our calculations at T = 
using the spin wave approximation. At finite tempera- 
ture, the long range order will be destroyed due to the 
Mermin-Wagner-Hohenberg theorem^, although short 
range correlations will persist. We now recall some prop- 
erties of the ground state of some unfrustrated Heisen- 
berg antiferromagnets, shown in earlier studies. 

1. Role of the dimension and the spin 

When the spins are quantum variables, the energy of 
the ground state is lower than the classical value, due 
to quantum fluctuations or "spin waves" in the case of 
periodic structures. As studies using a variety of analyt- 
ical and numerical methods have shown (see review^), 
the effects of quantum fluctuations vary depending on 
the dimension, on the spin quantum number S, and the 
type of structure. Quantum fluctuations are expected to 
become smaller as S increases, and as the dimension in- 
creases. It can be shown that, for hypercubic lattices, 
the spin wave expansion in fact is an expansion in 1/zS, 
where z = 2d. For d = 2 and S = h, as in the square 
lattice or the honeycomb lattice, it was an open question 
for some time whether quantum fluctuations would be 
strong enough to destroy long range order, until enough 
evidence^ 3 - was presented in favor of long range order. 

2. Role of the coordination number 

In a given dimension, and for simple lattices, one can 
ask what the local quantum fluctuations are when the 
coordination number is changed. In two dimensions, two 
examples of unfrustrated systems are the square lattice 
with z — 4 and the honeycomb lattice with z = 3 (see the 
review inM). Spin wave calculations, as well as Quantum 
Monte Carlo calculation o 13 ' 15 have shown that the order 
parameter is smaller in the z = 3 case. The values of the 



staggered magnetization, defined by the value of m s — 
\(S?)\ as obtained by QMC calculations are m s s q = 0.3173 
and rrig — 0.2788 (with sq and he standing for square 
lattice and honeycomb lattice respectively). This is in 
accord with the already remarked tendency towards a 
less classical behavior for systems of lower dimension, 
hence fewer nearest neighbors. 

In contrast, as pointed out ini£, when the coordination 
number is not constant, the local magnetization tends to 
be in fact larger when z is smaller. In the dice lattice, 
where sites can have z — 3 and z — 6, spin wave theory 
and QMC calculations have shown that it is the small z 
sites that have the larger value of the local staggered mag- 
netization, with m z s = 6 = 0.3754 while mf =3 = 0.4381. 
This result shows the "counter-intuitive" trend towards 
a less classical behavior for sites of bigger z in structures 
having a distribution of z values. Other systems showing 
this tendency include the quasiperiodic octagonal tiling 
and, as shown inii, the Penrose tiling. 



B. The Penrose tiling. 

We now describe the Penrose tiling, along with some 
of its properties. As shown in Fig[TJ which shows a finite 
portion of a perfect tiling, the tiles of this structure are 
a pair of rhombuses of angles 7r/5 (thin rhombus) and 
2ir/5 (thick rhombus). The edge length will be taken to 
be a = 1. 



1. Symmetries 

One remarkable feature is that despite the absence of 
perfect translational invariance, patterns of arbitrarily 
large size in different regions of the tiling can be made to 
overlap ("local isomorphism"). The mean repetition dis- 
tance of a pattern of linear size R is proportional to R, a 
property that replaces the strict translational invariance 
of crystalline structures. Rotational invariance holds in 
the same "weak" sense - for any given pattern, its equiv- 
alent under rotation by a multiple of 27r/5 can be found 
elsewhere on the tiling. 

The Penrose tiling possesses a hierarchical symmetry, 
being invariant under so-called inflation and deflation 
transformations. Inflation is a reversible operation which 
can be thought of as a set of decimations of vertices of 
the tiling, followed by a re-connection of the new ver- 
tices. The new tiling is defined on a length scale that is 
bigger by the factor r = v ^ +1 (the golden mean). The 
equivalence of the old and new tilings means here that 
no environments are created or destroyed in the process 
of inflation or deflation, and that one can find an exact 
match between any arbitrary (finite) regions of the two 
tilings. 

A rhombus-based structure such as the one shown in 
Fig[T] is a two-sublattice system. On the infinite tiling, 
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the two sublattices are equivalent, as regards the frequen- 
cies of the different vertices (see subsection below). 

2. Local environments 

As Fig[2] shows, many kinds of vertices are present in 
the tiling. In this paper, we have chosen to classify prop- 
erties of vertices according to their coordination number 
z. These site coordination numbers range from 3 to 7 in 
the Penrose tiling, with the average value z being exactly 
4. The figure shows the seven local environments possi- 
ble for coordination numbers z = 7,6,5,4 and 3. The 
first figure shows a five-fold symmetric site, which in fact 
comes in two varieties, the F (for football cluster) and the 
S (for star cluster). The properties of these sites under 
decimation are different. More will be said on these sites 
when we discuss the results for the ground state stag- 
gered magnetization distribution. We note that in the 
infinite tiling, each of the sublattices A and B have the 
same distribution of vertices of each type. 




FIG. 2: Local environments in the Penrose tiling. 



III. THE LINEAR SPIN- WAVE 
APPROXIMATION 

This unfrustrated Heisenberg antiferromagnet is ex- 
pected to have long range order at zero temperature, as 
remarked earlier. In the ground state, spins acquire a 
non-zero magnetization along the z-direction, with the 
total staggered magnetization of N = J2i( e i^i) 
where = ±1 depends on whether the site is on the 
A or B sublattice. 

Assuming that quantum fluctuations are small, one can 
expand around the classical antiferromagnetic ground 
state in which (eiSf) — S for all sites. We introduce, 
according to standard definitions (see for example^), 
Holstein-Primakoff boson operators for the sites on the 
two sublattices: ai,a\ for sites on the A sublattice, 
bi,b\ for sites on the B sublattice. The boson operators 
account for deviations from the classical configuration, 
with the spin deviation operator on A-sublattice sites 
given by afa, = n, = S — Sf and on B sublattice sites 
b^bj =nj = S + Sj. After linearization, the Hamiltonian 
is 

H LSW = -JS(S+l)N b +JS (4 a i+ b j b ]+ a l b ]+ b i a i) 

<ij> 

where Nf, = 2N is the number of bonds, N is the number 
of sites. This quadratic Hamiltonian can be diagonalized 
by a generalized Bogoliubov transformation, as shown by 
White, Sparks and Ortenburger— . Note that the diago- 
nalization can be carried out in the real space basis, and 
that there is no need to introduce collective operators as 
in periodic systems. 



3. The Penrose tiling viewed as a five dimensional object 

One standard method*- to obtain the Penrose tiling is 
by projecting a subset of vertices of a five dimensional 
cubic lattice onto a plane of a certain irrational orien- 
tation that is determined by the requirement of five-fold 
symmetry. As described in the Appendix, vertices are 
projected only if they lie within a five-dimensional cylin- 
der that runs parallel to the "physical" xy plane. This 
approach to generating the quasiperiodic tiling allows one 
an alternative visualization of the Penrose tiling in "per- 
pendicular space" - i.e. the vertices are projected onto 
the orthogonal three dimensional subspace. 

When the Penrose tiling is thus mapped into perpen- 
dicular space, one sees (as illustrated in the Appendix 
in Fig |20| . that the vertices now lie within four flat 
pentagon-shaped regions. Each family of sites of Fig[5] 
projects into a distinct domain - in other words, the 
perpendicular space representation allows us to separate 
sites according to their coordination number. We will 
use this useful property of the Penrose tiling in order to 
represent the complex antiferromagnetic ground state in 
a simpler way in Sec.V. 



A. Outline of numerical calculation scheme. 

In the expressions that follow, we number the sites 
such that the sites on sublattice A come first (i = 
1, ...,N/2), followed by the sites on sublattice B (i = 
N/2 + 1,...,N). The quadratic term in Hlsw can 
then be written in compact form in terms of the vector 
X T = (ai,a 2 , ....,a N/2 ,b\, ....,b\ /2 ) as 

Hlsw = JSX^H 2 X (2) 
where Hi is a real symmetric matrix of four ^ x T blocks: 



Z A 


C 


C T 


Z B 



(3) 



Za and Zg are diagonal matrices, with [Z a)u — %ii the 
coordination numbers on sublattice A, with a similar def- 
inition on sublattice B. The connectivites of the sites are 
given in the off-diagonal blocks, C: CV, = 1 if i and j 
are nearest neighbors and Cij — otherwise. Details of 
the numerical diagonalization using this real space ba- 
sis are described in£, where the LSW solution for the 
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ground state of the octagonal tiling was obtained. The 
main steps of the calculation will be outlined here for 
completeness. 

One seeks the canonical transformation taking the set 
{at, bj } to a set of boson operators {a m , f3 m }, in which 
the Hamiltonian is diagonal: 

N/2 

HlSW = Eq + ^ ^ti a ln a m + ^ m PLPm (4) 
m=l 

where Eq is the ground state energy 

N/2 

E = -JS(S + l)N b + JS n ™ (5) 

Two of the eigenmodes, apj/2 and 6jv/2i correspond to 
zero eigenmodes, and arise due to a global rotation of all 
the spins. These should be excluded from the diagonal- 
ization procedure, which transforms the remaining (N-2) 
operators. The generalized Bogoliubov transformation 
from the set {a, 6} to the set {a, (3} is X m = T mi X'i 
where X' T = (ai, a N/2 -i, p[,...., P N /2-i)- Tne 
rectangular T-matrix has the following structure: 

/ Ai,i .... A hN \ 



f-1,1 ■ ■ ■ ■ A f-l.N 

Bi,i .... Bi, N 



\ -Sf-ia ■ ■ ■ ■ Bn_ 1n J 

The transpose of the row vectors of this matrix, denoted 
by A m and B m , correspond to solutions of the eigenvalue 
equation 

gH 2 A m = il+A m 

gH2B m = —Q m B m (7) 

where the first equation holds for matrix elements of the 
rows in the upper half of T (the "positive subspace" ) , 
and the second for the rows in the lower half of T (the 
"negative subspace"). g is the matrix of commutators 

9ij = [x u x}], 

9ij = ±5ij (8) 

where the sign depends on whether i corresponds to the 
positive or the negative subspace. The new operators X' 
obey the same type of commutation relations [X'i, Xj] — 
g'ij — ±5ij (where i,j — 1, N/2 — 1) . 
The T-matrix should satisfy 

T^gT = g' (9) 

We note that in degenerate subspaces, a generalization 
of the Gram-Schmidt orthogonalization must be carried 
out with respect to the matrix g, in order to satisfy EqJ5] 



B. Finite size samples for numerical calculations. 

Our numerical solution of the linearized spin wave 
Hamiltonian is carried out on finite samples, using 
LAPACK 20 . We have mentioned in the previous section 
that quasiperiodic tilings can be considered as projec- 
tions of a simpler higher dimensional periodic systems. 
For our calculations, we have considered two types of 
samples. The first type are finite samples obtained by 
projection from 5D, for which we assumed open boundary 
condition. Such samples have the usual problem of spu- 
rious boundary modes. When considering averages over 
the entire spectrum, however, as one does in the compu- 
tation of the order parameter, we found good agreement 
between these finite samples, and periodic approximants 
described below. More precisely, spins deep in the inte- 
rior of these samples behave in the same way as the spins 
in the boundary-free periodic approximants, described 
next. 

In order to apply periodic boundary conditions, we 
have considered the so-called Taylor approximants of the 
Penrose tilingi&. These rectangular samples are obtained 
by slightly tilting the projection (or applying a "shear"), 
so that the resulting projection has a finite periodic 
length in the x- and y-directions, as shown in Fig|2T1 
These approximants clearly break the hierarchical and 
five-fold rotational symmetry of the Penrose tiling at 
length scales comparable to the repetition length. How- 
ever, one can generate Taylor approximants of bigger and 
bigger size, and thus approach the infinite Penrose tiling. 
Our samples have the additional property of having an 
equal number of sites in the A and B sublattices, thus 
ensuring a total spin of zero in the ground state. We 
considered approximants of 94, 246, 644, 1686, 4414 and 
11556 sites. 



IV. MAGNON SPECTRUM AND 
WAVEFUNCTIONS 

A. Energy spectrum and density of states. Linear 
dispersion law. 

The numerical Bogoliubov transformation gives the en- 
ergies as well as the T-matrix for each of the systems 
considered. The two sets of energies fi+ and become 
identical in the limit of infinite size when sublattices A 
and B become strictly equivalent. In our case, since the 
two sublattices are not identical, the numerical values 
are slightly different (less than a percent for the larger 
systems and not visible on the scale of our figures). 

Fig[3^) shows the integrated density of states (IDOS). 
We used the following expression for the IDOS, 



2 N/2 

N W ^nY, 6 ^- ^rn) (10) 
ni—1 
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N(E) 




0.2 0.4 0.6 0.3 1 1.2 1.4 



FIG. 3: a) (Color online) Integrated density of states N(E) 
versus E (expressed in units of J) calculated for three con- 
secutive approximants ( black dots : N = 644, red dots : 
N = 1686, grey : N = 4414) . b) Low energy tail of N(E) 
versus E 2 (the straight line is a fit to the data). 



The energy E is expressed in units of J, and the factor 2 
arises because the sum includes only half of the complete 
spectrum. Data for for three successive Taylor approx- 
imants has been plotted, with the points in light gray 
correspond to the largest system size, points in red to 
the medium size, and points in black to the smallest size. 
The integrated DOS is only slightly size dependent, as 
seen by the overlapping of the three sets of data points. 
The salient features of the figure are the several groups 
of closely spaced energy levels, the main gaps, which are 
stable with increasing system size, and a discrete jump 
at the energy E = 3, corresponding to a macroscopic 
degeneracy, as discussed below. 

The locations of gaps are more easily seen in the graph 
of the smoothed derivative o£N(E) - the density of states 
(DOS), which is shown in FigHJ One sees the character- 
istic fluctuating form of the density of states typical of 
quasiperiodic systems. 

One can distinguish several groups of energies sepa- 
rated by gaps. The highest energy bands around the 
values E « 5.4 and E ~ 6.4, correspond to wavefunc- 
tions that are localized on the 6- and 7-fold sites, as will 
be shown shortly. Lying below these in energy are the 
states centered primarily on the 5-fold sites. The peak 
at the exact value of E = 3 corresponds to string-like 
states living on the z — 3 sites. The lowest energies cor- 
respond to states of relatively extended character, as we 
will discuss below. 

As shown in FiglSJ}) the low energy part of the IDOS 
can be fitted to a power law of the energy, N(E) ~ E 2 
reflecting a linear dispersion of the magnon modes in 
this region of the spectrum. Fitting to a form N(E) = 
E 2 /(8irc 2 ) gives a sound velocity on the Penrose tiling 
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FIG. 4: Density of states calculated for the Taylor 6 approx- 
imant (N=11556 sites). 



of c = 1.08 J. It is interesting to compare this result 
for c with the corresponding values on the square lat- 
tice and the octagonal tiling, both having in common 
with the Penrose tiling the same value of the classi- 
cal energy per spin E c i/N = —2JS 2 . This value is 
c sq = 2^2 J Sa s» 1.41J on the square lattice (for edge 
length a = 1 and S — 1/2). On the octagonal tiling, our 
estimated value is c oc t a ~ 1.3 J, in agreement with Milat 
and Wessel's IDOS data£. To resume, an acoustic-type 
dispersion relation is obeyed at long wavelengths, with a 
spin wave velocity which is smaller than in the octago- 
nal tiling, which is in turn smaller than the value on the 
square lattice. This is presumably due to the fact that the 
density of sites is largest in the Penrose systems, followed 
by the octagonal tiling, and finally the square lattice. 

B. Wavefunctions. 

We now discuss the magnon wavefunctions and their 
spatial characteristics for the different parts of the en- 
ergy spectrum. We first show that the coordination num- 
ber plays an important role in determining the extent to 
which a site participates in the wavefunction for a given 
eigenmode, ip( E \ This can be seen from FigsJS^,) through 
e) , which show the weight fractions as a function of the 
energy, for each of the five values of the coordination 
number z. The weight fractions were defined as 

/n=El^ ) ! 2 /E^f ) i 2 C 11 ) 

where the i^,!^, ... are the subsets of sites whose coor- 
dination numbers are z = 3,4.... Only sites of sublattice 
A are considered (a similar calculation for sublattice B 
gives the same results). 

The plots of the weight fractions show a number of in- 
teresting properties of the magnon modes. In particular, 
they show the preponderance of certain types of site in 
the different energy bands. Thus it is clear from FigfS] 
e) that the highest energies correspond to states having 
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FIG. 5: a) - e) Weight fractions f n (see text) as a function 
of E for z = n (n = 3, 7), as computed for the Taylor 5 
approximant (N=4414 sites). 



a large amplitude on sites of z = 7. The highest energy 
band has a width of about 0.16 and is centred around 
E = 6.4. FiglT^,) represents one such state in real space, 
with sites are shown with varying intensity depending on 
the square of the wavefunction amplitude. The darkest 
spots are those on sites with z = 7. 

The group of states next highest in energy are those 
involving z — 6 sites, and correspond to energies in the 
range 5.43 < E < 5.49. These states have a smaller dis- 
persion than the states in the topmost group of energies, 
reflecting the fact that the six-fold sites are significantly 
fewer in number than the seven-fold sites. 

The middle band of energies, 3 < E < 5 arises for 
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FIG. 6: The average value of ipi(E) 2 per site for each z plotted 
against E for states in the bottom of the spectrum. 



states involving sites of z — 5, which account for about 
30% of the sites. From FiglS]:) one sees that the four- 
fold sites are particularly important in a narrow range of 
energy within this band. 

As seen by the large step in N(E) at E = 3, a large 
number of degenerate states occur at this energy. This 
is due to wavefunctions that have their entire support 
on the z — 3 sites. These are string-like states forming 
closed loops - on length scales that range from a small 
ring (around the footballs, for example) to being as large 
as the system size. A linear combination of such degen- 
erate states is shown in FigiTb). 

The lowest energy states , for E < 1, are the closest 
to extended states. The wavefunction amplitudes depend 
less sensitively on the site coordination number. The fact 
that all sites participate is best seen from Fig|5J where 
we have plotted for each of the five z-values the average 
probability for a given site, f n /N n , (where N n is the num- 
ber of sites of the nth family) , as a function of the energy 
E. Clearly, for the low energy states, the probability am- 
plitude is nonzero for all the values of z. On the other 
hand, not all of the sites participate in a given wavefunc- 
tion, and for a given energy the wavefunction is mainly 
confined to a set of disconnected patches. The lower the 
energy, the larger the patches where the wavefunction is 
non-zero. FigsEh) and d) show the states corresponding 
to two energies. For E — 0.8887 there are many closely 
spaced small patches. As the energy gets smaller, the 
patches of non-zero amplitude get larger, along with the 
spacing between them, and at energy E = 0.1692, for 
example, the wave function has two large patches (which 
are of course repeated periodically, due to the boundary 
conditions). 

To resume, the dimension of the support of eigenstates 
decreases as a function of the energy, from two at small E, 
to one at E = 3, tending to zero at the highest energies. 



Perpendicular space representation of the 
wavefunctions 



In the preceding section we alluded to the fact that the 
perpendicular space representation of the Penrose tiling 
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FIG. 7: Intensity plots representing the probability ipi(E)' 2 for 
different energies E (a darker shade corresponds to a higher 
probability): a) E = 6.469 b) E = 3.0000 c) E = 0.8887 d) 
E = 0.1692. The solutions have been obtained for the Taylor 
5 approximant (N=4414 sites). 



FIG. 8: Intensity plots in perpendicular space representing 
the probability ipi(E) 2 for different energies E (a darker shade 
corresponds to a higher probability): a) E = 6.469 b) E = 
3.0000 c) E = 0.8887 d) E = 0.1692 



is a useful way to see the environment-dependence of 
spatially varying quantities. We now illustrate this in 
the case of the four wavefunctions shown in Figs [7] rep- 
resented in perpendicular space by Figs[H Each vertex 
of the Taylor approximant is mapped (see the Appendix 
for more details) onto a point {x±, y±, z±}. We show 
the projection in the plane z± = 2 of the Penrose tiling, 
with regions shaded according to the local value of the 
wavefunction. Specifically, the intensity of the spot at 
site i is proportional to the value of \ipi(E)\ 2 . Fig(5Ji) 
shows the perpendicular space projection for the wave- 
function corresponding to the energy E — 6.469. The 
spots of maximum intensity are in the region that cor- 
responds to z = 7 (as seen in Fig |2"0)) . FigJSp) shows, 
similarly, that the wavefunction for E = 3 is non-zero for 
the region corresponding to z = 3. The last two figures 
show wavefunctions that are delocalized in perpendicular 
space (ie, all sites are involved, regardless of the value of 
the coordination number). 



D. Participation ratio and multifractal analysis. 

We present the results for the inverse participation ra- 
tio (IPR), defined by 



P-\E) = 



Recall that as N is increased, the inverse participation 



ratio decreases as 1/N for truly extended states, tends 
to a constant for localized states, and has an interme- 
diate behavior scaling as for the so-called critical 
states. This quantity has been much studied in the case 
of tight-binding models for electrons in disordered sys- 
tems, in particular, close to or at the critical disorder for 
the metal-insulator transition (MIT)22. For electrons in 
the Penrose and Ammann-Kramer-Neri tilings, Grimm et 
al22 have found a distribution of values of (3 ranging be- 
tween 0.5 and 1. Turning now to our spin problem, Fig [9] 
shows the results for log(P _1 (-E)) versus E for the Taylor 
approximants for three different sizes. (Note: the values 
have been calculated separately for each sublattice, ac- 
cording to the corresponding sector Or 11 of the eigenvalue 
spectrum). The fluctuations in the IPR tend to be quite 
large from one energy to the next, (note that the figure is 
plotted on the logarithmic scale), however, the smoothed 
IPR is an increasing function of E over most of the en- 
ergy spectrum. There is a noticeable dip in some of the 
IPR values as the energy approaches the value E w 3, 
and at E — 3 the value of the IPR does not reflect the 
spatial extent of the eigenstates, because of the mixing 
of the macroscopically degenerate states at this energy. 

Fig[9] shows the most marked size dependence in the 
region of small energies, in accord with our earlier ob- 
servation that states relatively delocalized at low energy. 
At high energies, the states are more localized, and the 
size dependence is accordingly smaller. 

As a general remark, results for the IPR should be 
treated with precaution, in the case of quasiperiodic sys- 
tems, as compared with disordered ones. In the latter 
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FIG. 9: Inverse participation ratio (IPR) plotted on a log 
scale as a function of the energy for three successive Taylor 
approximants 



case, away from the MIT, the decay of the wavefunctions 
is typically exponential, and degenerate states are un- 
likely. In the quasicrystal, exactly degenerate states do 
occur at special values of the energy {E = 3 here), which 
leads to an IPR value equal to that of extended states, 
due to linear combinations of localized states. Our results 
for the IPR, taken together with the preceding analy- 
ses of the wavefunctions, indicate that wavefunctions can 
be considered to be two dimensional, and power law ex- 
tended in the lower end of the spectrum. The average 
value of {3 , found by fitting the form P^ 1 (E) oc N~P 
found in this region is (3 w 0.9±0.1. 

The multifractal character of wavefunctions for 
quasiperiodic Hamiltonians can be considered to be es- 
tablished in a number of one dimensional models^ and 
strongly indicated in a number of two dimensional sys- 
tems including the tight-binding model for electrons 23 
and the phonon problem 2 ^, and even in three dimensional 
models^. The scaling properties of the wavefunctions 
can be investigated by calculating the f(a) spectrum. 
As f(a) depends on system size, and our system sizes 
are rather small, an extrapolation to infinite size is not 
without risk. Nevertheless, we have investigated the mul- 
tifractal scaling properties in the low energy end of the 
spectrum, where the states are extended, and have an 
patchy structure, with peaks and valleys spaced farther 
and farther apart as energy decreases - see Figs. [7b) 
and d). We find evidence for a nontrivial distribution 
of exponents in the limit of large system sizes. FigflOl 
shows the results of a multifractal analysis of two repre- 
sentative low energy states for three system sizes. The 
shape is a priori size dependent, and state dependent as 
well. One can nevertheless distinguish a smooth func- 
tional form of the f(a) function in each case. The curves 
have a maximum at ev(0) ~ 4 with /[a(0)] = 2.0, which 
is the dimension of the support - also called similarity 
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FIG. 10: Some examples of /(a) plots for small energy eigen- 
states in three Taylor approximants. 



dimension - of the wavefunction of this system. The so- 
called information dimension for these wavefunctions is 
given by /[a(l)] = a\ (also known as D\). It corresponds 
to the intersection of the f(a) curve and a straight line 
of slope 1 and we find cx\ ~ 1.5. For comparison, for a 
three dimensional system at the metal-insulator transi- 
tion (MIT) the critical values have been found to be^ 
a(0) = 4 and a(l) = 2. 



V. RESULTS FOR GROUND STATE ENERGY 
AND STAGGERED MAGNETIZATIONS 



A. Ground state energy 

The ground state energy per site is found from Eq[5] 
for the periodic approximants of size N ranging from 96 
to 11556. The energies are shown plotted against N~ 3 / 2 
in the figure QT] This is the expected power law for the 
ground state finite size correction in two dimensional pe- 
riodic systems. We find that it is obeyed on the aver- 
age in the Taylor approximants, with deviations that get 
smaller as size increases. Our extrapolation to infinite 
size gives an asymptotic value E — —0.643(0) ± 0.0001. 
The sign of the deviations varies with system size, and 
we believe this may be due to the defects present in the 
Taylor approximants. Since each approximant was ob- 
tained by a different rational section in five dimensional 
space, the distribution of local environments differ from 
sample to sample. In2I it was argued that the ground 
state energy of the tiling may depend in a simple way on 
moments of the variable zS. The first moment is just the 
average value of zS = 4S 1 in all cases. However ((zS)^ 1 ) 
varies from sample to sample, and this possibly accounts 
for the sign and magnitude of the observed deviations 
from the straight line power law behavior. 
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FIG. 11: The scaling of ground state energy of the Taylor 
approximants with system size 
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FIG. 13: Averages and standard deviations of the local stag- 
gered magnetization as a function of z for N=4414 ( LSW 
theory (dashed line) QMC (continuous line)) 




FIG. 12: A portion of a Taylor approximant with vertices 
colored according to the value of the onsite magnetization 
(magnetization values : red (small) , green (intermediate) blue 
(highest)) 
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B. Staggered local magnetizations 

The absolute value of the local magnetizations in the 
ground state are given in linear spin wave theory by 

m s (i) = \(Si)\=S- l A «| 2 (*<N/2) 

fe>f-i 

m s (j) = \(S*)\=S- Yl \ B ik? (i>iV/2)(12) 
fe<f-i 

Fig l 1 21 represents how local magnetizations vary in space 
on a portion of the Penrose tiling. The color of the cir- 
cles around each vertex varies from red (small magneti- 
zation) to blue (high magnetization). The lowest values 
of staggered magnetization are found on a certain subset 
of z = 5 sites, as will be discussed further below. The 
largest values are found on the low coordination sites of 

2 = 3. 

FigHSlindicates how values of the local magnetizations 
vary with z value for system size N — 4414. The figure 
shows the average values and the standard deviations of 
mi(z) for each value of z, as obtained in LSW theory 
and from QMC , as reported inii. Finite Penrose sam- 



FIG. 14: Probability distribution of the local staggered mag- 
netizations for the five different values of z (N= 11556) 



pies gave the same results as the periodic approximants 
when the sites on the free boundary layer were excluded 
from the analysis. The figure shows that Quantum Monte 
Carlo gives a narrower spread of the magnetization. Lin- 
ear spin wave theory clearly overestimates the fluctua- 
tions from average behavior, giving a too high value for 
Z = 3 and too low values for higher values of z. We show 
as an example, the details of the distribution of values ob- 
tained for z — 3 in Fig ll4l There is a continuum of values, 
as expected for a quasiperiodic structure, but also some 
pronounced peaks, corresponding to specific local config- 
urations. We explain these features in terms of the next 
nearest neighbor configurations present on the tiling, as 
seen in the next section. The five-fold sites come also in 
three main categories, as shown in Table 1, and this will 
be taken up in more detail further below. 

The size dependence of the magnetizations can be seen 
in FigfTSl where we show the values of the staggered 
magnetization in six different sample sizes (TV = 96 to 
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FIG. 15: System size dependence of the average local stag- 
gered magnetization for each coordination number z (lines 
indicate fit to power law iV -1 / 2 ) 



N = 11556). We show the scaling of the average value 
of rrii for a given z calculated for each of the six approxi- 
mants in Fig ll5l The averages are seen to obey a scaling 
with N^ 1 / 2 as in periodic systems. 

We now turn to the explanation of the fluctuations ob- 
served in Fig ll3l of the average staggered magnetization 
with z. This behavior can be explained with the help of 
a simple star cluster model, discussed in Sec. VI. 



FIG. 16: (Color online) Fig ll2l represented in perpendicular 
space. Two planes are shown (corresponding to z± — 1, 2) 
(magnetization values : red (small) , green (intermediate) blue 
(highest)) 



VI. PREDICTIONS USING THE HEISENBERG 
STAR MODEL. 

We first consider the effects of variation of the num- 
ber of nearest neighbors using a simple Hciscnbcrg star 
model, as first outlined in ie . A central spin So is coupled 
to z neighbors Sj. The Hamiltonian H = 52,- = i ^SoSj, 
can be expanded in boson operators a (describing the 
center site spin fluctuations) and bj, as described in 1 -. 
The corresponding magnetizations are respectively too = 
S — Smo and rrij = S — 8m j where (for z > 1) 



C. Perpendicular space representation of the 
ground state staggered magnetization 

We show, in Fig[T51 two sets of projections onto per- 
pendicular space of the vertices of the Taylor approxi- 
mants. The corresponding onsite magnetizations, as cal- 
culated for the tiling in real physical space, are colored 
according to their values as previously, in Fig[T2] The fig- 
ures show that, as expected, the domain corresponding 
to each coordination number has a distinct color. The 
colors are not absolutely uniform since no two sites are 
identical. Self-similar patterns can be investigated us- 
ing this type of representation, however our system sizes 
are not large enough to enable a quantitative analysis 
of self-similarity in the ground state. FigfTBb shows the 
domain corresponding to the F sites, as a central star- 
shaped region, which has the lowest m si values. The S 
sites project into a different star-shaped domain shown 
in Fig. \Wb , and can be seen to have a bigger value of 
the staggered local magnetization. We mentioned that 
the Penrose tiling is invariant under the inflation trans- 
formation where edge lengths are expanded by a factor 
r. If there is any self similarity of the ground state un- 
der an inflation operation, it can be perceived on the 
perpendicular space magnetization map as pairs of simi- 
lar patterns in regions which are related via an inflation 
transformation. 



Smo — l/( 2 — 1) 
Smj = l/z(z — 1) 



(13) 



Quantum fluctuations of the center spin are thus larger 
than those of the outer spins. This is a consequence of 
the fact that in such clusters, the classical term, which 
creates an onsite potential Vo oc z, thereby discourag- 
ing boson formation on the center site, is dominated by 
the transverse terms. These terms in a)b^ create and 
annihilate boson pairs on each link, leading to quantum 
fluctuations being greater for sites with more neighbors. 

This cluster model is a first approximation. In order to 
explain the multiple peaks in the values of m s ; seen nu- 
merically for a given z, we must take into account longer 
range structural details. 

We consider therefore the two tier Heisenberg star 
shown in FigfTTl as described inii. The Hamiltonian 
of this cluster is linearized after introducing Holstein- 
Primakoff operators a^,ai, (i = 1, ...,zz') and bj, (j = 
1, .., z) The resulting expression for the center site mag- 
netization is 



m s (z, z) 



1 



fi{z,z>)-zft{z,z>)-Az>- 



(14) 



where / 1(2) = -z' ± (2 - z + y/A - Az + (z + z') 2 ). 

The main new feature of m s (z,z') is its non- 
monotonicity. to is shown plotted against z for various 
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FIG. 17: Two-tier Heisenberg star 



values of z' in Fig[TH] For each of the curves of fixed z', 
a shallow minimum is seen to occur for values z ~ 1 + z' , 
that is, when the coordination numbers of the center site 
and the nearest neighbors are matched. 

Returning to the Penrose tiling and approximants 
thereof, we use the number of next nearest neighbor 
bonds n(z) to classify sites of a given z. The z = 3 sites 
can be classified in three subgroups, depending on z', as 
shown in Table 1. Fig [HI shows the values of m s (z,z') 
obtained for the Penrose tiling sites. To take the case 
of z = 3, three values corresponding to three main local 
configurations are found, to be compared with the three 
principal peaks of the red curve in FigO The substruc- 
tures arise due to differences in the third nearest and 
further neighbor configurations. 

The number of next nearest neighbors also serves to 
distinguish between the different z — 5 sites, of which 
there are three main types. There are two varieties of 
five-fold symmetric sites: the football-shaped clusters (F) 
and the star-shaped clusters (S). The former have z' = 2, 
while the latter have z' — 4. The F sites thus have 
the smallest onsite magnetizations and the S-sites, on 
the contrary, have the largest onsite magnetization. The 
remaining (most frequently occurring) z = 5 sites which 
do not have a local five-fold symmetry have intermediate 
values of z' , 



z 


z' (corresponding frequency) 


m s {z, z') 


3 


4 (31%),4.33 (27%), 4.67 (42%) 


0.41, 0.42, 0.43 


4 


3 (100%) 


0.36 


5 


2 (14%),2.4 - 3.2 (81%), 4 (5%) 


0.26, 0.35, 0.41 


G 


3 (100%) 


0.37 


7 


2.3 (100%) 


0.33 



Table 1. Values of z (and their frequencies) for each 
coordination z and the predicted values of m(z,z') using 
EqE] 



A. Fourier transformation 

Fourier transforms of the original structure, as well as 
that of the antiferromagnetic ground states are shown in 



m(z, z r } 



0.45 




3 4 5 6 7 

FIG. 18: m a (z,z') plotted as a function of z (see text) for 
selected values of z' (as given in Table 1). 



V 



-4 -2 2 4 

FIG. 19: Intensity plot of the static longitudinal magnetic 
structure factor £?"(k) for the N = 4414 Taylor approximant. 
The highest intensity peaks are shown in the figure at posi- 
tions that are indicated by circles(squares) for the magnetic 
(nonmagnetic) structures respectively. The relative intensity 
is denoted by a linear gray scale ranging between zero (white) 
and maximum intensity (black). 



1191 The highest intensity peaks are shown in the figure 
at positions that are indicated by circles (squares) for the 
magnetic (nonmagnetic) structures respectively. Peaks 
can be indexed using four indices which are integers for 
non-magnetic peaks and integer or half-integer for the 
magnetic peaks. As in the case of the octagonal tiling 6 
the half integer indices arise due to the doubling in size of 
the antiferromagnetic unit cell in five dimensions. This, 
coupled with the extinctions in the structure factor rule 
due to the multiplicity of each unit cell, then leads^ to 
the observed "shifting" or "displacement" of the mag- 
netic peaks with respect to the nonmagnetic ones. 



VII. CONCLUSIONS 

We have presented a detailed analysis of energies and 
wavefunctions of magnon eigenmodes found by numeri- 
cal diagonalization of the Heisenberg model on approx- 
imants of the Penrose tiling in the linearized spin wave 
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approximation. A linear dispersion is found at low ener- 
gies, as in other bipartite antifcrromagnets, despite the 
lack of translational invariance. This is to be expected 
in the long wavelength limit, where structural details do 
not play a role, as known in related studies of low energy 
vibrational modes in quasicrystals or in glasses. The ef- 
fective spin wave velocity as obtained from the low en- 
ergy spectrum was found and compared with that in the 
square lattice and the octagonal tiling. 

The spatial dependence of eigenmodes has been inves- 
tigated. It is shown that in different ranges of energies, 
states have their support on sites of a different given co- 
ordination number. At low energies, in contrast, all sites 
participate, and the resulting eigenstates are relatively 
extended, as seen by the absolute values and the size de- 
pendence of the participation ratio. A few low energy 
extended eigenstates are analysed for their multifractal 
scaling properties, including the similarity dimension and 
information dimensions. In general, our studies indicate 
that the dimensionality of the states diminishes progres- 
sively as energy increases. 

We find the ground state energy of this antiferromag- 
net, and give the distribution of the local magnetizations 
as a function of site coordination number. Perpendicular 
space projections are shown to demonstrate the simplic- 
ity of the ground state in this representation. The struc- 
ture factor of the magnetic state is found and compared 
with that of the nonmagnetic state. 

A simple analytical model, the two-tier Heisenberg 
star, is presented to explain our results. The role of next 
neighbors is shown to be nontrivial, and explains the z 
dependence of the onsite staggered magnetizations. 
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where 9 — 2ir/5. 

For a vertex to be selected, its perpendicular coordi- 
nates must fall within the rhombic icosahedron W which 
is the projection of the 5D unit cube. In view of the ex- 
pression for M± given above, W decomposes into four 
subdomains W i} corresponding to z± = 1,2,3, or 4 
mod(5). The acceptance window Wi in each of these four 
planes in the perpendicular space is a pentagon. Two 
such pentagons are shown in Fig |20l (corresponding to 
z± = 1 and 2), the other two being the same upto an 
inversion. 

The parity of the vertex as given by ni determines 
the sublattice to which it belongs. Thus, the points that 
project into the planes z± = 1,3 correspond, say, to sub- 
lattice A , while the planeszj^ = 2,4 correspond to sub- 
lattice B. In the infinite tiling, the two sublattices are 
equivalent, and the "even" and "odd" windows are the 
same upto an inversion. 

Finite samples of arbitrarily large size can be readily 
generated using this selection criterion, by considering 
large enough volumes of the 5D cubic lattices. 

Turning to the question of the local environments, we 
show in Figf20] a portion of the Penrose tiling, with ver- 
tices of different z colored differently. When these ver- 
tices are represented in perpendicular space, the colored 
domains are clearly seen corresponding to each value of 
z. We have seen that there are three types of z = 5 sites, 
including the F and S subsets. The F sites project onto 
the central pentagonal regions in the planes z± = 1 and 
4, while the S sites project into the central pentagons in 
the planes zj_ = 2 and 3. 



APPENDIX A: OBTAINING THE PENROSE 
TILING AND ITS APPROXIMANTS BY 
PROJECTION 

1. The perfect Penrose tiling. 

The perfect Penrose tiling can be obtained by the pro- 
jection onto the xy plane of selected vertices of a five 
dimensional (5D) cubic lattice. A vertex, designated by 
the five-dimensional vector n±, n 2 , n5, is selected for 
projection if its projection in the three dimensional "per- 
pendicular space" belongs in the region W (called the 
"selection window"). The projection matrices Mu and 
M± give the real space coordinates x\\ , y\\ and perpendic- 
ular space coordinates x± ,y±, z±_ respectively are 



2. The Taylor approximants. 

To obtain periodic approximants, the physical plane 
has to be tilted so that it has a rational orientation (that 
can be chosen arbitrarily close to the original irrational 
orientation). There will then be two lattice vectors Ay 
and A 2 of the 5D lattice, whose projections give the pe- 
riodic lengths in the physical plane. This is achieved 
in practice by using an oblique projection M± appro- 
priately defined depending on the choice made for A\ 
and A 2 - The resulting selection window W will be a de- 
formed rhombic icosahedron. Details of the construction 
of the projection operators for the so-called Taylor ap- 
proximant can be found ir>i£. This method can be used 
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(a) (b) 




FIG. 20: (Color online) (upper) Portion of the tiling showing 
vertices colored differently, according to coordination number 
z (lower) The same tiling after projection into perpendicular 
space (a) the plane z± = 1, b) the plane z± = 2 ) 




FIG. 21: Taylor r 3 approximant (N=644 sites) 
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